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Abstract This study used realistic representations of cloudy atmospheres to assess errors 
in solar flux estimates associated with ID radiative transfer models. A scene construction 
algorithm, developed for the EarthCARE mission, was applied to CloudSat, CALIPSO and 
MODIS satellite data thus producing 3D cloudy atmospheres measuring 61 km wide 
by 14,000 km long at 1 km grid-spacing. Broadband solar fluxes and radiances were then 
computed by a Monte Carlo photon transfer model run in both full 3D and ID independent 
column approximation modes. Results were averaged into 1,303 (50 km) 2 domains. For 
domains with total cloud fractions A c < 0.7 top-of-atmosphere (TOA) albedos tend to be 
largest for 3D transfer with differences increasing with solar zenith angle. Differences are 
largest for A c > 0.7 and characterized by small bias yet large random errors. Regardless of 
A c , differences between 3D and ID transfer rarely exceed ±30 W m 2 for net TOA and 
surface fluxes and ±10 W m~ 2 for atmospheric absorption. Horizontal fluxes through 
domain sides depend on A c with ~20% of cases exceeding ±30 W m~ 2 ; the largest values 
occur for A c > 0.7. Conversely, heating rate differences rarely exceed ±20%. As a cursory 
test of TOA radiative closure, fluxes produced by the 3D model were averaged up to 
(20 km) 2 and compared to values measured by CERES. While relatively little attention 
was paid to optical properties of ice crystals and surfaces, and aerosols were neglected 
entirely, ~30% of the differences between 3D model estimates and measurements fall 
within ±10 W m~ 2 ; this is the target agreement set for EarthCARE. This, coupled with the 
aforementioned comparison between 3D and ID transfer, leads to the recommendation that 
EarthCARE employ a 3D transport model when attempting TOA radiative closure. 


H. W. Barker (El) 

Cloud Physics and Severe Weather Research Section (ARMP), Environment Canada, 
4905 Dufferin St., Toronto, ON M3H 5T4, Canada 
e-mail: howard.barker@ec.gc.ca 

S. Kato 

NASA — Langley Research Center, Hampton, VA, USA 


T. Wehr 

European Space Agency, Noordwijk, The Netherlands 


<£) Springer 


658 


Surv Geophys (2012) 33:657-676 


Keywords Cloud ■ Radiation • Climate ■ Satellite ■ EarthCARE ■ CloudSat ■ Cloud- 
Aerosol Lidar and Infrared Pathfinder Satellite Observations ■ The Moderate Resolution 
Imaging Spectroradiometer • Clouds and Earth’s Radiant Energy System 


1 Introduction 

The ultimate boundary conditions of Earth’s climate system are defined by the flow of solar 
radiation in and longwave emission out of the top-of-atmosphere (TOA). All else consti- 
tuting climate and life falls between these boundaries. Understandably then, radiative 
transfer through, and optical properties of, the Earth-atmosphere system lie at the core of 
our concepts of climatic forcing and feedbacks. Moreover, many remote sensing tech- 
niques, both passive and active, involve measurement and subsequent model interpretation 
of (largely atmospheric) radiative transfer. Mathematical treatments of radiative transfer 
were developed initially by astrophysicists (e.g., Chandrasekhar 1950), but the latter half of 
the 20th century saw substantial advancements from Earth scientists in their response to the 
myriad conditions and scales of Earth’s climate system. Arguably, these advancements are 
exemplified best by work centred on radiative transfer for cloudy atmospheres. 

Over the past three decades there has been a string of diagnostic studies aimed at 
demonstrating errors associated with neglect of 3D radiative transfer within cloud, weather 
and climate models (e.g., McKee and Cox 1974; Davies 1978; Ellingson 1982; Welch and 
Wielicki 1984; Stephens 1988; Davis et al. 1990; Barker and Davies 1992; Cahalan et al. 
1994; Marshak et al. 1995; O’Hirok and Gautier 1998; Barker et al. 1999, 2003; Cole et al. 
2005a). Generally, these studies used small, and often selective, samples of 3D cloud fields 
derived either from idealized models, such as arrays of simple forms (Welch and Wielicki 
1985; Kobayashi 1988) and scaling algorithms (Marshak and Davis 2005a), passive 
satellite data (Cahalan et al. 2005), aircraft data (Barker 1992; Raisanen et al. 2004) or 
limited-area domains simulated by cloud system-resolving models (Barker et al. 1999; 
Hinkelman et al. 2005). Ideally, a large number of 3D cloud fields based on observational 
data should be used to demonstrate errors due to neglect of 3D solar transfer. This would 
enable proper hypothesis tests to be carried out involving clouds whose geometric prop- 
erties resemble most faithfully those that actually occur. 

Cloud properties derived from synergistic retrievals using data from lidars, cloud-pro- 
filing radars and passive radiances provide, thus far, the most complete pictures of true 
cloud structure (ESA 2001; Stephens et al. 2002; Hogan et al. 2006; Winker et al. 2007; 
Dupont et al. 2011). However, because current active-passive systems point fixedly at 
either zenith (ground-based) or nadir (satellite-based), they provide 2D vertical cross- 
sections, not 3D domains. Those derived from satellite data have the advantage of pro- 
viding a global perspective on a regular grid whose horizontal spacings are not subject to 
height-dependent advection rates. Their disadvantage is that they are less resolved than, 
and not as sensitive as, their surface counterparts. 

To get around the issue of satellite retrievals being essentially 2D cross-sections, 3D 
fields can be constructed using the algorithm of Barker et al. (2011) that was developed for 
use by the EarthCARE satellite mission (ESA 2001). EarthCARE is slated for launch in 
2015 and will carry a high-spectral resolution lidar, a 94-GHz cloud Doppler radar, a 
7-channel imager and a broadband radiometer. Their algorithm combines conventional 
passive satellite imagery with nearby coeval information from active-passive profiles; 
recipient columns associated with off-nadir passive pixels get assigned, by proxy, an 
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active-passive donor column associated with a nadir pixel thereby leading to production of 
a 3D domain of cloud. 

This study represents a first step towards comparing ID and 3D solutions for large 
samples of cloud fields derived from satellite data. A-train data and retrieval products 
(Kato et al. 2010) construction algorithm developed by Barker et al. (201 1) to construct 3D 
cloudy atmospheres. A broadband solar transport solver (capable of affecting ID and 3D 
solutions) was then applied to them. Reflected fluxes from the 3D model were compared to 
corresponding (20 km) 2 CERES values (Wielicki et al. 1996). For (50 km) 2 domains that 
resemble the size of weather and climate model grid-cells, flux differences arising from ID 
and 3D solutions were compared. 

In the following section the A-train dataset is described briefly. This is followed by 
descriptions of the 3D construction algorithm, the radiative transfer model and the nature 
of the simulations. Results are presented in Sect. 6, and conclusions and recommendations 
are in the final section. 


2 Merged A-train Data 

Data gathered by several A-train satellites were used to assess the importance of 3D solar 
radiative transfer. A merged A-train dataset has been compiled at NASA-Langley in which 
CloudSat (Stephens et al. 2002) and CALIPSO (Winker et al. 2007) profiles were mapped 
to 1-km resolution thereby associating each column with a 1-km MODIS pixel (Kato et al. 
2010). It also includes the line of (20 km) 2 CERES broadband radiances (and fluxes) along 
the CloudSat-CALIPSO track as well as MODIS pixels either side of it out to ~ 40 km. 

A merged cloud mask was created by interpolating CloudSat’ s and CALIPSO’ s cloud 
masks onto CERES’s arbitrary vertical grid which includes temperature, humidity and 
ozone (from the Goddard Earth Observing System — Data Assimilation System — GEOS-4; 
see Bloom et al. 2005; Kato et al. 2005). If a CloudSat or CALIPSO cloudy cell overlaps a 
CERES layer, the cell is designated to be filled with cloud. 

Because a large fraction of columns have just cloud masks from CloudSat and CALI- 
PSO, vertical profiles of cloud properties were set using MODIS quantities and the fol- 
lowing algorithm. For each MODIS pixel along the CloudSat-CALIPSO track, effective 
visible optical depth T eff , effective particle size r eff and particle phase were retrieved 
based on inversion of a ID radiative transfer model (Minnis et al. 2008, 2010a, b). If there 
are n ice and n liq cells in a column, the former having temperatures <273 K and the latter 
>273 K, it is assumed that 


where X and £ are ice and liquid water contents, respectively, p x 0.917 g cm~ 3 , and A Zi 
is layer thickness. If CloudSat’ s columnar classification (Sassen and Wang 2008) is cirrus, 
altostratus, altocumulus or deep convection, all ice cells in the column use r ice = r e ff, and 
any liquid cells use r| iq =10 pm. When the columnar classification is stratus, stratocu- 
mulus, cumulus or nimbostratus, all liquid cells use r liq = r eff , and any ice cells use 
r ice = 50 pm (e.g., Yang et al. 2007). Equation 1 is then solved for water contents 
assuming that X = £ but restricting X to be no larger than 0.3 g m~ 3 (e.g., Protat et al. 



(1) 


2007). 
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This simple algorithm produces clouds that are, for the most part, locally homogeneous 
in the vertical. Obviously this can influence both ID and 3D radiative transfer solutions, as 
well as their differences which are of concern here. Nevertheless, when clouds produced 
this way are operated on by ID shortwave and longwave radiative transfer algorithms, 
local TOA flux estimates agree well with CERES data, particularly longwave values 
(Barker et al. 201 1). 

Without defending this vertical allocation scheme too much, it is worth noting that for 
geometrically thin clouds that get resolved into a small number of layers, such as many 
boundary layer clouds, the impact of vertical structure on 1D-3D transfer differences is a 
minor issue. The same goes for very cloudy scenes, which as will be seen constitute the 
majority of the cloudy (50 km 2 ) domains. Finally, preliminary studies, not reported here, 
using data from cloud system-resolving models suggest that vertical homogenization of 
water for towering convective clouds generally exaggerates 1D-3D differences, implying 
that results presented below are approximately an upper-bound for 1D-3D differences. 


3 Constructing 3D Cloud Domains from A-train Data 


For the convenience of readers, Barker et al. (2011) construction algorithm is reiterated in 
this section. A-train data make-up the 2D retrieved cross-sections (RXS). Let all pixels/ 
columns along the RXS at positions (i, 0), where for this study i = 1 was in the southern 
hemisphere of an ascending portion of an orbit, form the set of potential donor columns 
that can act as proxies for off-nadir recipient columns associated with pixels located at 
{ (i, j) | / £ [—7, — 1] U [1,7]} . For each recipient pixel, begin by computing 


F(U j- = 


k= 1 L 


n{i, j) ~ r k (m , 0)' 


: m £ [i — m \ , i + m 2 ] 


( 2 ) 


where r k are MODIS radiances, and k denotes spectral interval. For this study 
mj = m 2 = 200 pixels which corresponds to searching for ~±200 km along the RXS 
beginning at the pixel on the RXS that is in the recipient’s cross-track line. Next, F(i, /; m) 
are ordered from smallest to largest to form the set 


{min[F(i, j; m )], . . ., max[7 7 (i, j\ m)]} = {F(i,j\ 1), . . .,T(i,j\nn + m 2 + 1)}: 

F(i,j;n)^J r (i,j:n+ 1). 


( 3 ) 


Defining Euclidean distance between a potential donor at (m,0) and recipient at ( i , j) as 

D(i, j ; m) = A L\J (i - m) 2 +j 2 , (4) 

where AL is imager resolution, let D(i, j;m) and m go passively along with the ordering of 
F(i, /; m) so that {7 r (i. /; n)} has an associated co-ordered set of distances {V(i,j; n)}. One 
then solves for the index tn as 


argmin : f £ (0, 1), (5) 

m* e [1 , (mi +m2 + 1 )/] 

which means: find the index m that corresponds to the smallest distance between recipient 
and those pixels that constitute the smallest 100 f% of F(i, j; m) . For this study / = 0.03 so 
that the smallest 3% of the (/«[ + m 2 + 1) values of F{i,j\ m) are considered (provided 
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they have the same surface type as the recipient). Knowing m leads directly to m and thus 
the donor column. Finally, all the properties associated with the column at ( m , 0) get 
replicated at ( i , j). In addition to cloud properties, this includes profiles of temperature, 
moisture and aerosol, as well as surface conditions. The vast majority of cross-sections 
used here were completely over ocean and so there were 133 layers from the surface up to 
65 km. 

Figure 1 shows a schematic of this process which gets applied until the desired 3D 
scene is constructed. Clearly, RXS columns (j = 0) identify themselves so for them there 
is no need to apply the algorithm; the RXS forms the centre of the constructed domain. 

Verification of the construction algorithm using A-train data is not straightforward. One 
can, however, go a certain distance by attempting to reconstruct the RXS itself. In so doing, 
one attempts to fill an RXS column at (;, 0) by searching the RXS and applying (5) to 
potential donor pixels in [;' — m\, i — n] U [; + n, i + m 2 ] which bars the first ±n pixels 
next to i ; hence defining a dead-zone in the search process. This test is meant to mimic 
filling off-RXS columns that are ±n pixels away from the RXS. For example, when n = 5, 
searching for a proxy column begins five pixels away, just as for off-RXS pixels at (i, ±5). 

Figure 2 shows attempts to reconstruct a 400-km-long stretch of RXS. The upper image 
is the actual RXS merged cloud mask. This is an especially demanding case as it involves 
fairly dense multi-layer clouds. Over much of this domain passive-only retrievals would 
yield very little, if any, useful information about cloud vertical structure. Lower images are 
reconstructions for discrete values of n = 1, 5, 10 and 20. These results stem from using 
four spectral channels (0.62-0.67; 2.105-2.155; 8.4-8. 7; and 11.77-12.27 pm) in (2). By 
n = 5, which corresponds to the outer edge of an 1 1-km-wide domain the likes of those to 
be computed for EarthCARE, it is clear that some error was creeping in; nevertheless, a 
significant amount of detail was captured. Even out at n = 20, multi-layers of clouds were 
replicated well. The region where the greatest difficulty was encountered between 100 and 
200 km along the horizontal where clouds were transitioning or dying out entirely. 

Figure 3 shows mean profiles of several variables accumulated out to n as functions of n 
for the field shown in Fig. 2. Results for accumulated fields, of widths 2 n + 1, are shown 
because the algorithm is intended to produce full 3D domains not single rows. For these 
accumulations averaging included the original RXS as it is included in constructed fields, 
yet gave double weight to the reconstructed lines so as to represent scene construction on 
both sides of the RXS. For clouds higher than 10 km, layer cloud fraction A c , mean cloud 



Fig. 1 The thin RXS is shown along with the sequence of MODIS visible pixels associated with it. The 
objective is to fill the volume marked by the wider MODIS swath with cloud properties drawn from the 
RXS. For example, the column associated with the pixel at (in, 0) has been designated as the proxy for the 
pixel at (/, j) so the cloud-radiation attributes associated with (m, 0) get donated to (/'. j). The algorithm is 
applied until all desired off-RXS pixels are filled by donor RXS columns 
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Fig. 2 Topmost image is 
merged cloud mask (1-km 
horizontal resolution) for a 
stretch of tropical RXS along an 
orbit on 19 April 2007. Black and 
grey indicate ice and liquid, 
respectively. Sequences of lower 
images represent corresponding 
masks produced by the 
construction algorithm for 
various dead-zone lengths as 
listed 


4 channels 
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water contents (w), and v = ((' w)/a w ) 2 , where a w is standard deviation of w, were 
reproduced extremely well for all n. The largest errors on A c occurred for n = 20 at 
altitudes near 5 km. (w) contents were reconstructed very well except for n = 5 and 10 for 
clouds between 2 and 5 km high where associated values of v were overly small, indicating 
that these errors were the result of having selected a small number of columns with 
anomalously large w. On the other hand, for n = 20, almost all reconstructed clouds below 
10 km lacked sufficient horizontal variability as indicated by values of v being twice as 
large as they should be, despite corresponding (w) being fine. This was due to too many 
occurrences of a small number of RXS columns being used multiple times. 



Fig. 3 Domain-average profiles of layer cloud fraction, cloud water content (w), and v = ((w) /^J 2 , where 
c t w is standard deviation of w, for the RXS (actual) shown in Fig. 2 as well as for reconstructed cross- 
sections corresponding to various dead-zone lengths as listed 
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As mentioned, the exemplary field used here was a difficult case of multi-layer tropical 
cloud. Most 400-km sections of cloud do not exhibit as much intricacy as this one. 
Numerous other examples were examined and almost all reconstructions performed 
equally well or better than those shown here. In general, the more extensive and planar the 
clouds and the fewer the number of definite layers, the better the reconstruction. 


4 Radiative Transfer Model 

Broadband shortwave (SW) fluxes and radiances were computed by a 3D Monte Carlo 
algorithm (Barker et al. 1998, 2003). For this study, full 3D results pertain to horizontal 
grid-spacings AL ~ 1 km. Cyclic horizontal boundary conditions were employed. Inde- 
pendent Column Approximation (ICA) results were produced by the same model using 
AL = 10 8 km. Gaseous transmittances (H 2 0, C0 2 , O 3 ) were computed using the correlated 
L-distribution method with 3 1 quadrature points in cumulative probability space (Scinocca 
et al. 2008). Optical properties for liquid droplets (Dobbie et al. 1999; Lindner and Li 
2000) and ice crystals (Fu 1996; Fu et al. 1998) were resolved into four bands. Aerosol 
effects were not included but Rayleigh scattering was. For each grid-cell a cumulative 
distribution of extinction was computed including Rayleigh scattering, gases and cloud 
particle sizes segregated into 23 radius bins progressing in width following a power law 
from 0-1 pm to 128.0-161.3 pm. Size-integrated Mie scattering functions were used for 
each bin, and phase functions for ice crystals were treated as though they were spherical 
droplets. To reduce variance in radiance estimates, Mie functions were used for a photon 
packet’s first 5 scattering events with the Henyey-Greenstein function, with proper 
asymmetry parameter, thereafter if needed (Barker et al. 2003). 

The majority of data used here were collected over the Pacific and Indian Oceans. 
Surface reflectance for oceans was represented by the Cox and Munk (1956) ergodic wave- 
facet model. Probability of reflection and direction of reflected photons were derived from 
zenith angle of incident photon and surface wind speed as reported in the A-train dataset 
(Bloom et al. 2005). 

Land was treated as Lambertian with spectral albedos set to values inferred from 
MODIS data as reported in the A-train dataset (Rutan et al. 2009). For the model’s 
0.2-0. 7 pm band, values obtained from 0.469 pm radiances were used. For the three near- 
IR bands, values for 1.24 pm were used. 


5 Radiative Transfer Simulations 

Radiative transfer simulations were performed on constructed domains that measured ~61 
km wide by ~ 14,000 km long and consisted of 85.4 x 10 6 columns and ~ 113.6 x 10 6 
cells. Going much wider with the current version of the construction algorithm is not 
recommended (see Figs. 2, 3 as well as Barker et al. 2011). Less memory was required 
than one might think as optical properties were unique only for RXS profiles and each off- 
nadir profile was indexed to its RXS donor. While photons showered entire domains 
uniformly, the outer 5 km in the across-track direction and 500 km at the along-track ends 
were omitted from the analysis. The omitted portions did, however, serve in the 3D 
simulations as buffer-zones that at once shielded the inner-domain from potentially adverse 
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effects set-up by cyclic horizontal boundary conditions and facilitated horizontal transport 
of photons through the sides of the inner-domain. 

Simulations were performed using solar geometry at satellite time-of-crossing. As such, 
solar azimuth angle relative to A-train ground-track <p 0 (i, j ) and solar zenith angle 9 0 (i, 
j) were set to values reported in the merged A-train dataset. Photons were injected uni- 
formly across constructed domains and given an initial weight of j) = cos 9 0 (i, j). 


6 Results 

Two lines of investigation were explored and these are discussed here. First, TOA reflected 
fluxes inferred from CERES radiances were compared to their simulated counterparts 
which were obtained by applying the 3D radiative transfer model to constructed A-train 
domains. Second, TOA reflected solar fluxes obtained by running the radiative transfer 
model in 3D mode were compared to results obtained by running it as a ID model. 

6.1 Reflected TOA Fluxes: CERES Versus Modelled 

Comparison of CERES quantities to modelled values constitutes an attempt to perform a 
radiative closure experiment as CERES data were not used in the retrieval process. Fig- 
ure 4 shows the CloudSat quick-look image for granule 999. While the constructed 3D 
domain was 14,000 km long and 61 km wide, the dots on the track run from ~44°S to 
~72°N and demarcate the portion of the orbit considered for comparison. 

Figure 5 shows TOA fluxes inferred from CERES radiances (Loeb et al. 2005, 2006) as 
well as model estimates for the constructed domain. Each value corresponds to a (20 km 2 ) 
CERES pixel. 256 x 10 6 photons were injected over the full domain making for 
~ 120 x 10 3 per CERES footprint with maximum Monte Carlo noise errors for TOA 
reflectance of ~ 0.002. Visually, the agreement is strikingly good for the most part, and 
encouraging too given the lack of synergistic retrieval and relatively little attention paid to 
ice scattering properties and surface albedo and no attention given to aerosols. 

Figure 6 shows the cumulative histogram of differences between modelled-CERES 
fluxes (neglecting values between 50°N and 60°N where some cloud-related data were 
missing). Only about 30% of the differences fall within the ±10 W m~ 2 goal of Earth- 
CARE (ESA 2001; Domenech et al. 2011) and the bias for this particular orbit is —15 W 
m~ 2 . This large bias, however, arises primarily from the near-cloudless, high-Sun cases 
between 5°S and 25°N. Underestimation of albedo in this region could result from 
undetected small clouds (that were not included in the model), systematically too small 
ocean albedo (stemming from underestimated surface wind speeds) or neglect of aerosol. 
The plot also shows the cumulative histogram of differences between modelled-CERES 
nadir radiances scaled by n in order to be of similar magnitude to fluxes. The fact that 
histograms for fluxes and radiances are almost identical suggests that, for this orbit, there is 
no clear advantage to using either nadir radiances or fluxes to assess the retrievals. It is 
worth pointing out, however, that while radiances are a direct measurement CERES fluxes 
carry the additional error of radiance-to-flux conversion. It is expected that inclusion of off- 
nadir radiances, as in EarthCARE, will make for a more challenging assessment of 
retrievals. 

Note also that measured radiances, and hence ID MODIS retrievals, are subject to 
radiative smoothing and roughening (e.g., Marshak and Davis 2005b). Thus, when they are 
used at the 1 km scale to initialize a radiative transfer model, even a 3D model as here, 
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July 5 2006 (CloudSat granule 999) 
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Fig. 4 Upper panel shows a CloudSat quick-look image for one of the orbits used here. The line indicates 
the satellite’s ground-track and the dots on it mark the start (~44°S) to the end ( ~ 72°N) of the portion used 
for radiation calculations. Solar geometry is shown in the lower left while solar zenith 0 o and azimuth <p 0 
angles along with surface wind speed for the marked area on the quick-look image are shown on the graph 


secondary smoothing and roughening takes place. This undoubtedly introduces errors into 
the comparison that are far from straightforward to estimate. The hope is that column 
properties retrieved by a robust synergistic algorithm will rely less on passive retrievals 
and thereby reduce smoothing/roughening errors. There is also the possibility that an 
expanded rendition of the construction algorithm might partially address the effects of 
horizontal transport (cf.. Barker et al. 2011). 

Before leaving this section it is interesting to give an example that goes well beyond the 
initial steps of radiative closure using domain averages and consider the much more 
demanding case of imagery reconstruction. Figure 7 shows a 50 x 150 km domain of 
MODIS cloud optical depth t and its construction algorithm counterpart. Next to this is the 
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TOA reflected SW flux 



Fig. 5 CERES and 3D Monte Carlo estimates of TOA reflected solar fluxes for the marked portion of the 
orbit shown in Fig. 4. The 60-km-wide MODIS 0.645 pm image is shown along the top 


Fig. 6 Cumulative frequency 
distributions of (3D model— 
CERES) TOA reflected solar 
fluxes for the data shown in 
Fig. 5. Also shown is the 
distribution of differences 
between CERES nadir radiance 
and 3D Monte Carlo nadir 
radiances (multiplied, for plotting 
purposes, by n so as to put them 
on roughly the same scale as 
fluxes) 


July 5, 2006 



corresponding 0.645 pm MODIS image and nadir reflectances produced by the 3D model 
using constructed i. While the reconstructed radiances are not perfect, the main features 
were captured well and will only be that much better once full synergistic retrieval comes 
online. 

6.2 Reflected TOA Flux Differences: ID Versus 3D 

So, while a serious attempt at radiative closure using the current A-train dataset would be 
premature, on account of known difficiencies in retrieved properties it is likely that cloud 
structural features are defined well enough to conduct a comparison between ID and 3D 
mean solar fluxes for mesoscale-size domains. In so far as the RXSs and their 3D con- 
structed counterparts capture realistic cloud structure down to the 1 km scale with 
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Fig. 7 Image on the extreme left (actual) shows cloud visible optical depths inferred by applying an inverse 
ID RT model to MODIS radiances for a 50 x 150 km scene. Next to it is its 3D constructed version. To 
their right are MODIS radiances at 0.645 pm, and next to it corresponding radiances from the 3D RT model 
using the constructed field. Solar zenith angle was 40° and the Sun was coming in from the direction shown 


sufficient accuracy, results presented here constitute a fair assessment of radiometric errors 
due to neglect of 3D solar transport. 

When comparing domain-average characteristics of ID and 3D radiative transfer it is 
common to employ standalone domains that are assuredly isolated in space via cyclic 
horizontal boundary conditions. This is because one knows from the start that all photons 
either get absorbed by the domain’s atmosphere or underlying surface, or exit through the 
top. An arbitrarily delineated domain in the real atmosphere, however, is not isolated, 
cyclic boundary conditions do not apply and a clean comparison between ID and 3D 
transfer becomes complicated — the more so as the domain’s horizontal extent diminishes. 
For this reason, particular attention was paid in these experiments to fluxes through the 
vertical sides of domains. 

Results for each simulation were averaged into 50 x 50 km sub-domains, a nominal 
area representing weather and climate model grid-cells. As each simulation used 
256 x 10 6 photons, each sub-domain received ~750 x 10 3 photons. This translates into 
very small Monte Carlo errors for domain-average fluxes and heating rates, ~2.5 times 
smaller than those for results presented in the previous sub-section. 

Results were partitioned into sub-domains with total cloud fractions A c \ (i) 
0.05 < A c < 0.3; (ii) 0.3 < A c < 0.7; and (iii) A c > 0.7. This is because previous afore- 
mentioned studies have suggested that cloud structural effects should maximize at inter- 
mediate cloud amounts and taper off for larger and smaller A c . Data from five orbits were 
used and the number of 50 x 50 km sub-domains for the three ranges of A c were: (i) 205; 
(ii) 246; and (iii) 852, for a total of 1,303 sub-domains. 

Figure 8 shows flux differences (3D- ID) as functions for /t 0 for the three ranges of A c . 
For small cloud fractions the only clear signal is greater reflected flux at TOA when 3D 
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transfer was considered. The same trend is visible for A c e [0.3,0.7]. This is the well- 
documented effect of cloud-side illumination for /to <0.8 (Welch and Wielicki 1985). 
Conversely, at large /t 0 there is a slight tendency for 3D transfer to produce smaller 
reflectances due to photons emerging from cloud-sides in predominantly downwelling 
directions. Additionally for A c e [03,0.7] there are trends for larger atmospheric absorption 
and weaker surface absorption (irradiance) by 3D transfer as )( 0 decreases. Again, this is due 
to increased cross-sectional area of cloud presented to direct-beam as )t 0 decreases. The vast 
majority of differences for TOA reflected and surface absorption are smaller than ±20 W 
m 2 ; for atmospheric absorption, differences are largely confined to ±5 W m 2 . 

When cloud fractions become greater than 0.7, the dominant feature is a large scatter of 
positive and negative differences that cluster around zero with a very slight tendency for 
3D transfer to yield smaller reflectances by about 3 W m~ 2 and larger atmospheric 
absorption by ~2 W m 2 . These values are not too dissimilar to those obtained by Cole 
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Fig. 8 Differences between fluxes predicted by 3D and ID transfer for (50 km 2 ) scenes constructed from 
A-train data. Left column is for reflected flux at TOA, centre is for atmospheric absorption, while the 
rightmost is for surface absorption. Each row corresponds to scenes with A c as listed on the plots. Grey lines 
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Fig. 9 Cumulative frequency 
distributions for differences in 
reflected solar flux at TOA as 
predicted by a 3D transport 
model and its ID counterpart for 
50 x 50 km constructed scenes 
(i.e. the values shown on the left 
column of Fig. 8) 


reflected flux at TOA 



et al. (2005b) using 2D model atmospheres, and possibly stem from photons undergoing 
internal reflections between multiple cloud layers thereby increasing photon pathlengths, 
number of scattering events, and hence absorption. As f 0 decreases, the direct beam can 
undercut layers of clouds when 3D transfer is admitted; for ID transfer, however, photons 
must diffuse through upper clouds before undergoing internal reflections. 

Figure 9 shows cumulative frequency distributions of differences in reflected flux at 
TOA predicted by 3D and ID radiative transfer for all sub-domains partitioned into the 
three cloud fraction ranges. Median values for A c < 0.7 are very close to zero but for 
cloudier scenes it is —8 W m 2 . It can be expected that errors in TOA reflectance will 
exceed ±10 W m~ 2 for 10-30% of scenes having A c < 0.7 if ID rather than 3D transfer is 
used. This rate jumps to ~ 50% of scenes with A c > 0.7. Preliminary results suggest that these 
percentages increase by 15-20% when smaller domains the size of CERES footprints are 
considered. The implication therefore for EarthCARE is that 30-65% of its constructed 
scenes can expect errors in computed reflected short wave (SW) flux at TOA to exceed 
its ±10 W m~ 2 target if ID rather than 3D transfer is adhered to. More dramatic errors 
in excess of ±50 W m~ 2 appear as though they will be very rare regardless of cloud 
conditions. 

Figure 10 shows surface and atmospheric absorption histograms that correspond to 
those in Fig. 9. Regarding surface absorption, the most interesting point is that errors due 
to use of ID transfer are often large (i.e. greater than 20 W m 2 ) for A c > 0.3 and are 
dominated by overestimation relative to 3D values. This stems from no illumination of 
cloud-sides in ID calculations, especially apparent for A c e [0.3, 0.7]. Biases in total 
atmospheric absorption due to use of ID models appear to be of secondary importance for 
(50 km 2 ) domains as they rarely exceed ±10 W m~ 2 . Naturally more extreme localized 
differences exist at smaller scales but they are largely averaged out at 50 km. 


6.3 Horizontal Transport Through Domain Sides 

For a column of atmosphere associated with a 50 x 50 km sub-domain, let R , T and A 
be net fluxes out the top, at the surface and absorbed in the column. Since the 


<£) Springer 


670 


Surv Geophys (2012) 33:657-676 


surface absorption 


0.9 

0.8 

0.7 


■A <0.3 
c 

■ 0.3 < A <0.7 
c 

0.7 < A 


// 


flux difference (W m ) 


atmospheric absorption 



1 1 1 1 1 1 1 1 1 

1 1 1 1 




A 

c 

0.3 

0.7 


/ s * 

/ s* 



< A <0.7 

/ 



<A 

,7 





/ 

/ 






( 











1 

1 






1 

1 






/ 






s'J 






1 1 1 1 

i i i i 

i i i i 



-20 -10 0 10 20 
flux difference (W m" 2 ) 


Fig. 10 As in Fig. 9 except these distributions are for atmospheric and surface absorption (i.e. data shown 
in the centre and rightmost columns of Fig. 8) 
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Fig. 11 Schematic diagram illustrating horizontal transport H of solar radiation. The column of concern is 
in the centre of each figure and consists of an overcast cloud. R, A and T denote reflectance, atmospheric and 
surface absorptance, respectively. Left figure shows ID results where there is no transfer through the 
column’s vertical sides and so H iD = 0. The other figures show the impact on H 3D due to clouds outside the 
column of concern 


sub-domains considered here were effectively standalone and not subject to cyclic 
boundary conditions, 

R T + A 7 ^ ti 0 F 0 ( 6 ) 

holds in general where Fq is incident flux on a normal surface at the TOA. Horizontal 
transport of radiation through the sides of a sub-domain can be defined as 

H = l i 0 F 0 -{R + T+A). (7) 

Marshak and Davis (2005b) provided an in-depth discussion of H for solar transfer through 
cloudy atmospheres with (7) corresponding to their (12.13). Figure 11 illustrates how to 
interpret (7). For ID transfer ( 6 ) is an equality and so H = 0. For 3D transfer, however, if 
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Fig. 12 As in Fig. 9 except 
these distributions are for implied 
horizontal transport of radiation 
as defined in (7) 



the shadow of a cloud outside the sub-domain is cast into the sub-domain, but the majority 
of reflected photons come from, for example, a large planar cloud aloft (so that Ri D ~ 
Rid), H > 0 which might seem odd given that transport through the column’s sides is less 
than it would have been had the external cloud been absent. If the external cloud is now 
moved so that direct-beam radiation undercuts the layer cloud and some radiation gets 
scattered by the external cloud back into the sub-domain, H < 0. Since non-zero H can 
arise through differences between ID and 3D values of R, T or A as well as actual 
horizontal flux through vertical boundaries, perhaps it is best to think of H as implied 
horizontal flux. 

Figure 12 shows histograms of H for the three classes of sub-domains using data from 
the five orbits considered here. Values of j/f| >50 W m~ 2 occur for ~ 10% of the scenes 
with A c > 0.7 and point directly to the potential importance of 3D effects on the radiation 
budget of atmospheric volumes resembling those of global model grid-cells. Although the 
scale at which H was evaluated here (50 x 50 km domains) and the nature of the domains 
differ much from those assessed by Marshak and Davis (2005b), values are similar to 
those they arrive at (cf. their Fig. 12.4). When sub-domain size is reduced to that of a 
CERES footprint, distributions of H for cases with A c e [0. 3,0.7] and A c > 0.7 resemble 
one another closely but the fractions of occurrences of |7/|>50 W m~ 2 increase to 
roughly 20%. Again then, it appears to be crucial that 3D solar transfer be utilized by the 
EarthCARE mission if it is to make a serious go at both estimating surface-atmosphere 
radiation budgets and assessing the quality of active-passive retrievals by demanding that 
simulated TOA fluxes be within ±10 W m~ 2 of values inferred from its broadband 
radiometer measurements for (10 km 2 ) scenes (Domenech et al. 2011), especially for 
broken cloud scenarios. 

Finally, Fig. 13 shows sub-domain average heating rates computed by the 3D model 
for a 13,000 km stretch of the orbit on 5 July 2006. It also shows percentage differences 
between 3D and ID transfer. The thing to point out is that differences rarely 
exceed ±20%. When they do, however, the bias persists from the surface to the base of 
those clouds that are responsible for maximum heating. Once a larger number of orbits 
get processed the intention is to examine heating rates for scenes sorted according 
to various cloud properties including cloud-type (e.g., Xu et al. 2005; Bacmeister and 
Stephens 2010). 
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Fig. 13 Beneath the 60-km-wide visible MODIS image are solar heating rate profiles computed by the 3D 
Monte Carlo model using 50 x 50 km sections of the 13,000-km-long constructed domain for the orbit 
shown in Fig. 4. The lower panel shows percentage differences between solar heating rates predicted by the 
3D and ID transport models 


7 Conclusions and Recommendations 

The primary purpose of this study was to begin to get a global picture of errors in estimates 
of solar fluxes, for cloudy domains roughly the size of those found in global models, that 
are incurred by applying the ID independent column approximation (ICA). Both the ICA and 
3D benchmark values were obtained by applying a Monte Carlo photon transport model to 
cloud fields constructed from A-train satellite data. The results are of relevance to both global 
models that use ICA-based methods to compute mean radiative fluxes (e.g., Randall et al. 
2003; Morcrette et al. 2008; Shonk and Hogan 2008) and to radiative closure studies. 

3D cloud scenes were created from 2D A-train satellite data via application of the 
Barker et al. (2011) 3D scene construction algorithm. Their algorithm was developed for 
the EarthCARE mission which intends to perform continuous radiative closure experi- 
ments to assess the quality of its active-passive retrievals. As many A-train cloud profiles 
consisted of a cloud mask only, profiles of cloud properties had to be allocated subject to 
constraint by MODIS retrievals. Based on preliminary results using cloud system- 
resolving model (CSRM) data, it appears that the simple technique used here slightly 
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overestimates differences between ID and 3D transfer models. Hence, one way view of 
the results is that they represent an approximate upper-bound for differences between ID 
and 3D solutions. 

Other limitations of this study were its: inability to address impacts of horizontal 
fluctuations in cloud extinction at scales smaller than ~ 1 km; treatment of scattering by 
ice crystals as though they were liquid spheres; and neglect of aerosols. It did, however, 
include important fluctuations at scales larger than 1 km (e.g., Zuidema and Evans 1998), 
and according to simulations by Cole et al. (2005b), cloud albedo depends weakly on 
horizontal resolution for grid-spacings less than ~ 2 km. The general impression, there- 
fore, was that 3D cloud structural features, plus some information pertaining to micro- 
physical variations, were inferred well enough from A-train data to begin setting the stage 
for a global statistical analysis of differences between ID and 3D solar transfer. The 
quality of this comparison is expected to improve once greater attention is given to ice 
cloud, aerosol and surface optical properties, and synergistic active-passive retrieval 
methods mature (e.g., Hogan et al. 2006). 

A Monte Carlo transport algorithm was applied to constructed scenes using A-train 
time-of-crossing values of cosine of solar zenith angle f.i 0 . Calculations were performed on 
domains measuring 61 x 14,000 km and results were partitioned into 1303 50 x 50 km 
sub-domains which nominally represent the size of cells found in conventional global 
models. In corroboration with previous studies, 3D-1D differences are largest for sub- 
domains with intermediate values of total cloud fractions A c , 3D transfer leading to larger 
reflectances at TOA and atmospheric absorption, and smaller surface irradiances. Results 
were, however, somewhat dependent on fi Q with mean differences being almost zero at 
large /.i 0 and about ~ 10 W m~ 2 for TOA reflectance, — 10 W m~ 2 for surface irradiance, 
and ~5 W m~ 2 for atmospheric absorption for )i 0 ~ 0.2. Similarly, as shown for a 
13,000 km stretch of a single orbit, CERES-inferred and model-predicted solar TOA fluxes 
agreed to with ±10 W m~ 2 about 30% of the time. Similar values were found for other 
orbits. 

More so than biases the most defining characteristic of 3D-1D flux differences, 
regardless of both fi Q and A c , was large random variations. For instance, for A c > 0.3 it was 
not unusual to find surface irradiance differences exceeding ±30 W m~ 2 . This scatter of 
results even occurs for near-overcast conditions. The reason for this may be related to 
multiple layers of clouds that lead to larger photon pathlengths, and thus slightly larger 
atmospheric absorption, but can sway reflectance and surface irradiance deviations one 
way or another. 

It is difficult to say whether deviations of these magnitudes will have significant ram- 
ifications on simulations of clouds, weather and climate. The few experiments with CSRMs 
that have attempted to explore deviations like these suggest not (Mechem et al. 2008; Cole 
et al. 2005a; Pincus and Stevens 2009). But those tests were far from exhaustive. More 
telling, perhaps, are the Monte Carlo Independent Column Approximation (McICA) 
experiments with large-scale models (Pincus et al. 2003; Barker et al. 2008; Morcrette 
et al. 2008). In those studies the amount of radiative noise associated with random sam- 
pling of clouds far exceeded that shown in Fig. 8. Nevertheless, there was very little impact 
on the evolution of simulations. 

What these results do suggest is that in order to limit errors in radiative closure studies, 
regardless of whether TOA fluxes or radiances get used, one should at least conditionally 
apply a 3D solar radiative transfer algorithm to retrieved/constructed scenes that exhibit 
highly variable cloud. This is especially so if one seeks to assess retrieval performance for 
small individual scenes using solar measurements. This is indeed the plan for the 
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EarthCARE mission which, as far as the authors are aware, will represent the first attempt 
in the atmospheric sciences to employ a 3D radiative transfer model operationally. 
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